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I.  INTRODUCTION 

We  are  interested  in  determining  validated  sets  of  elementary 
chemical  reactions  for  use  in  predictive  combustion  models.  For  this 
purpose,  a  simulation  of  the  laminar,  premixed,  one  dimensional,  steady 
state  flame  has  been  Implemented.  Iliis  has  the  advantage  of  being  a 
relatively  simple  combustion  problem  that  also  yields  predicted  tempera¬ 
ture  and  species  profiles  that  can  be  compared  with  suitable  burner 
experiments. 

The  basic  idealism  is  of  an  infinite  column  of  premixed  gas.  At 
some  point  the  gas  is  ignited.  A  flame  front  forms  and  propagates  along 
the  column  of  gas.  All  effects  are  assumed  to  be  one-dimensional. 
Eventually,  the  flame  will  reach  steady  state.  A  basic  characteristic 
of  the  solution  is  the  flame  speed,  the  velocity  with  which  the  flame 
front  moves  with  reference  to  the  undisturbed  gas  mixture.  The  flame 
front  is  marked  by  steep  teiiq>erature  and  chemical  species  profiles.  The 
above  formulation  is  referred  to  as  an  unbounded,  or  adiabatic  flame. 

In  practice,  a  flame  is  usually  stabilized  on  a  burner.  That  is, 
the  flame  will  propagate  to  the  burner  surface.  Loss  of  heat  to  the 
burner  will  stabilize  the  flame.  This  creates  different  boundary 
conditions. 


This  report  describes  a  method  of  numerically  solving  the  equations 
governing  a  laminar,  premixed,  one-dimensional,  steady  state  flame. 

Both  unbounded  flames  and  burner  stabilized  flames  are  discussed.  To 
solve  such  systems,  we  have  modified  a  standard  package  for  integrating 
one-dimensional  partial  differential  equations  so  that  it  efficiently 
handles  flame  equations.  The  basic  procedure  is  to  integrate  the 
equations  in  time  until  the  steady  state  solution  is  reached.  Changes 
have  been  made  in  the  basic  structure  of  the  code,  as  well  as  adding  a 
number  of  subroutines.  The  program  has  been  implemented  for  the  ozone 
flame  and  the  H2-O2-N2  flame. 

Section  II  describes  the  partial  differential  equations  governing 
the  flames.  In  Section  III,  the  basic  code  used  to  solve  the  equations 
is  introduced.  Section  IV  describes  the  modifications  made  in  the  code, 
and  Section  V  the  additional  subroutines  added  to  the  program.  Section 
VI  discusses  the  changes  in  the  equations  and  the  code  necessary  to 
model  burner  stabilized  flames.  Section  VII  reports  some  of  the  numerical 
considerations  in  successfully  using  the  code.  Finally,  the  ^>pendix 
gives  a  complete  listing  of  the  code,  plus  the  output  for  a  run  modeling 
a  typical  H2-O2-N2  flame. 
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II.  THE  UNBOUNDED  FLAME  EQUATIONS 


The  equations  for  a  multicomponent  reacting  ideal  gas  mixture  can 
be  found  in  the  literature. We  are  interested  in  the  equations  that 
describe  a  one-dimensional,  laminar,  premixed  flame  that  propagates  in 
an  unbounded  medium.  The  effects  of  radiation,  viscosity  and  body  forces 
are  ignored.  The  momentum  equation 


p 


h. 

A 

at 


3x 


.  3£ 

A 

ax 


can  be  eliminated.  Here  v  is  the  fluid  velocity,  p  is  the  density,  and 
p  is  the  pressure.  The  independent  variables  are  time  t  (sec)  and  dis¬ 
tance  X  (cm).  For  flames^j  the  fluid  velocity  is  much  less  than  the _ 

speed  of  sound,  that  is,  v^  «  p/p.  For  steady  state  flow,  this  implies 
that  the  pressure  gradient  is  negligibly  small.  So  the  pressure  is 
assumed  to  be  constant. 

The  pertinent  equations  are  then,  overall  continuity. 


at  ax 


(1) 


That  is,  any  change  in  density  with  respect  to  time  is  due  to  the  overall 
convection  of  the  mixture,  which  equals  the  gradient  of  the  total  mass 
flux  pv. 


Continuity  of  species: 

A  AA  *  AA  a  a 

p  =  .  pv  — i  ^  (pY  V.  )  +  R^M^,  k  =  1,2,...N 

at  ax  ax 


(2) 


^F.  A.  ffittiamBs  Combustion  Theorut  Addison-Wsley ^  Reading,  MA,Chapt9r 
1.  1965. 

2 

R.  B.  Bird,  W.  E.  Steuart  and  E.  S.  Lightfoot,  Transport  Phenomena, 
John  Wiley  and  Sons,  NY,  Chapter  18,  1960. 

R.  M.  Fristram  and  A.  A.  WesterQjerg,  Flame  Structure,  MoCfraa-Hill,NY, 
Chapter  V-1,  1965 

^J.O.  Hirahf elder,  C.  F.  Curtiss  and  R.  B.  Bird,  Molecular  Theory  of 
Gases  and  Ljauida,  2nd  Printing  corrected,  with  notes,  John  Wiley  emd 
Sons,  NY,  Chapter  11.1,  1964. 
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where  Yj^  Is  the  mass  fraction  of  the  k  species,  Vj^  is  its  diffusion 

velocity is  its  molecular  weight,  and  R.  is  the  rate  at  which  the 
species  is  produced  or  consumed  by  the  chemical  reactions.  So  changes 
in  the  concentration  of  the  species  can  be  due  to  the  convection  of 
the  species,  the  diffusion  of  the  species  in  the  mixture,  9;  {he  product¬ 
ion  or  consumption  of  the  species  by  chemistry.  The  term  pY^Vif  is  the 
mass  flux  of  tte  k^h  species^ (relative^to  v).  Note  that  the  total  mass 
flux  for  the  k^”  species  is  p(v  Vj^)  Yj^. 

Conservation  of  Energy. 


-P 


N 

1 

k-l 


V  Vk 


(3) 


A  A  At. 

where  T  is  the  temperature,  c  .  is  the  specific  heat  of  the  k  species, 

A  pK  ^ 

c  is  the  specific  heat  of  the  mixture,  X  is  the  thermal  conductivity  of 

*  A  ♦K 

the  mixture,  and  h]^  is  the  specific  enthalpy  of  the  k  species.  So 
changes  in  the  temperature  can  be  due  to  the  convection  of  heat,  the 
conduction  of  heat,  the  production  or  consumption  of  energy  by  the 
chemical  reactions,  and  a  small  amount  due  to  the  diffusion  of  species 
with  different  specific  heats.  Hie  have  not  written  a  negligibly  small 
term  due  to  the  species  gradients  (Dufour  effect) . 

The  thermal  equation  of  state  is  given  by  the  ideal  gas  law 


p  ■  P  R  T 


N 

2: 

k=l 


(4) 


where  R  is  the  gas  constant.  The  caloric  equation  of  state  is 


(S) 


th 


where  h.  °  is^the  specific  enthalpy  of  the  k'""  species  at  some  reference 
temperature  t  The  specific  heat  of  the  mixture  is  given  by 

O  e 


i 


7 


f  1  f 


^  •  I  f 

■!  I 


c  «  ,  c  ,  Y, 
p  k«l  pk  k. 


Froa  conservation  of  mass  we  have  the  relation 


i  ^k  “  1 

k»l  * 


N  *  “ 

.2,  ‘k\  =  0. 

k»l 


The  boundary  conditions  are  the  following.  For  x  ■  - 


T  .  Ty  and  (k-1,2, . . .N) . 


where  is  the  temperature  and  the  Yj^^  are  the  mass  fractions  of  the 
original,  undisturbed  mixture.  For  x  «  •* 

T  -  Tg  and  Y^  •  Y^^,  k-l,2,...N  .  CIO) 

A  A 

where  Tg  is  the  temperature  and  are  the  mass  fractions  of  the  burned 

mixture.  T.  is  called  the  adiabatic  temperature.  Since  we  are  assuming 
no  heat  loss  to  the  surroundings,  Tg  and  Yj^g  depend  only  on  the  chemistry, 
and  can  be  calculated  in  advance. 

In  practice,  the  Integration  will  be  over  a  finite  interval 

(x, ,  x_),  where  the  flame  front  will  be  located  roughly  in  the  center  of 
L  K  A 

the  interval.  The  boundary  conditions  (9)  will  be  applied  to  x, . 

However,  the  conditions  (10)  are  not  convenient.  After  the  flame  front, 

there  is  a  long  recombination  period  before  the  mixture  reaches  adiabatic 

conditions.  Choosing  the  interval  of  integration  long  enough  to  cover 

the  entire  recombination  zone  would  be  con^utationally  expensive.  So 

normally  the  flame  will  be  cut  off  before  it  reaches  adiabatic  conditions. 
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and  the  weaker  boundary  conditions 


37  3k  n  u  1  1  M 

-s-  “  — A-  ■  0,  k»l,2,...N  , 

3x  3x  Cll) 

A 

will  be  applied  at  x^^  . 

A 

At  this  stage  the  partial  differential  Eq.  CiJ  for  the  density  p 
can  be  eliminated  by  introducing  a  new  coordinate  $  such  that 

A 

5?  (x.t)  =  p(x"  ,  t)dx'  .  (12) 

J  Y 


_  3!1»  * 

Then  -i-  =  p 

and 

3t|i 

A 

3x 

3t 

notation  Eqs. 

(2) 

and 

a* 

!!il  = . 

m 

o 

3“^ 

3;i; 

^  A  A  A 


A  A  A 


'X*0 


'\A/  ^ 


^  (PVk)  *  \  Vp 


and 


'V* 

8^  °  a?  2  9j!V  ) 


N 

Z 

k-1 


3^ 


T  ^  3T  , 

Pk  k  k  a 

k»l  ^ 


(13) 


(14) 


^  ^  y  2 

where  the  tilde  variables  are  functions  of  t  »  t  (sec)  and  v  (gm/cm  ). 

For  numerical  convenience,  dimensionless  forms  of  Eqs.  (13)  and  (14) 

are  integrated.  That  is,  we  will  define  T  ■  'i>/To,,  t  «»  ^/t»  and  i)i  « 
where  t.  and  4).  are  chosen  so  as  to  obtain  reasonable  time  and  space 
scales  for  a  given  flame,  and  T.  is  chosen  so  that  T  is  the  same  order 
of  magnitude  as  the  larger  Yj^'s  .  Then  the  equations  are 


3Y, 


3Y.. 


t„  3 


—  ^  \  ‘  CIS) 


and 
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m  II  *  i*  r  1  3  1  1 

♦;:  “o  w  *  ;;;7  3^  ^ 


(16) 


N  N 

I  R^  y(»  T.  )  -  I  CJ,^  Y^V^  |I  )  , 

k=l  k=l 


where  the  variables  are  functions  of  t  and  t|'.  The  boundary  conditions 
are 


.  Y^  =  Y^u  •  1^  =  1.2...N 


at  <|i  =  4^^  and 

II  ,  !!ii 

3^  3i|) 


k  a  1.2. ..N 


(17) 


(18) 


at  T|»  * 

The  integration  can  start  from  any  initial  profiles,  as  long  as 
there  is  enough  energy  available  to  begin  the  combustion.  As  shown  in 
Section  V,  the  mass  flux  m  through  the  origin  is  modified  as  the 
integration  proceeds  so  as^to  keep  the  flame  front  in  the  center  of  the 
interval  of  integration  (4>, ,  i|)„) .  The  integration  proceeds  until  m 

L  K  .y  O 

approaches  a  constant,  and  the  terms  K  .  3T  approach  zero.  The 

3t  3t 

equation  for  species  N  is  not  integrated;  Y^  is  found  from  the  relation 
(7). 


The  burning  velocity,  S  ,  can  be  calculated  at  steady  state.  From 
the  continuity  equation  3(pvy  /3i|»*0,  or  pv  is  constant  with  respect  to 
\|i.  Then  we  can  take  any  of  the  equations  (15)  and  integrate  over  any 
interval  (a,b)  to  obtain 


10 


pvj  [Y^(b)  -  Y^Ca)]  =  / 


W 


P  (-»)  fY^(b) 


Vkl  ^ 


At  steady  state,  all  of  the  v.  are  equal.  Since  in  the  present  coordinate 
system  the  flame  does  not  move,  v.  (-»)  is  the  speed  at  which  the 
unburned  gas  is  approaching  the  flame.  Conversely,  S  =  Vj^  (-“)  is  the 
speed  at  which  the  flame  is  propagating  into  the  unburned  mixture. 

The  parameters  of  the  equations  (15)  and  (16)  are  the  specific 
heats  Cpj^,  the  specific  enthalpies  hj^,  the  chemistry  production  terms 

Rj^,  and  the  diffusion  velocities  Vj^  for  each  species,  plus  the  thermal 

conductivity  X  of  the  mixture.  The  c^j^  and  hj^  are  functions  only  of  the 

temperature,  and  can  be  evaluated  very  accurately  using  sixth  degree 
polynomial  fits^*  The  chemistry  is  generally  the  least  well  known  of 
the  input  data.  We  need  to  know  which  species  are  involved,  which 
reactions  can  occur,  and  the  rate  constant  for  each  reaction.  The  rate 

constants  k.  will  be  of  the  form  a  ^  ^  exp  (c/^) ,  where  either  b  or  c 
can  be  zero^  We  can  then  find  the  rate  r.  for  each  reaction  by  multiply¬ 
ing  the  rate  constant  times  the  concentrations  of  the  reactants  (con¬ 
centrations  of  the  k^h  species  =  p  Y^^/  Mj^).  Each  Rj^  is  then  found  by 

adding  the  rates  of  the  reactions  in  which  the  species  is  a  reactant. 
The  transport  parameters  X  and  Vjj  are  in  general  very  complicated  func¬ 
tions  of  temperature  and  species  concentrations.  Because  of  the 
numerical  complexity,  various  approximations  to  these  quantities  have 
been  used.  We  have  the  capability  to  run  the  code  with  a  number  of 
different  levels  of  approximation.  This  is  discussed  in  detail  in 
another  paper®. 


^S.  Gordon  and  B.  J.  McBride^  "Computer  'Program  for  Calculation  of 
Complex  Chemical  Equilibrium  Compoeitione ^  Rocket  Performance, 
Incident  and  Reflected  Shocks,  cmd  Chapman- Jouguet  Detonations", 
NASA-SP-273,  1971,  (1976  program  version) . 

^T.  P.  Coffee  and  J.  M.  Heimerl,  "Transport  Algorithms  for  Premixed, 
Laminar,  Steady-State  Flames",  to  be  published  in  Combustion  and 
Flame. 
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For  future  reference,  we  will  discuss  here  the  simplest  transport 
algorithm  used  in  the  code.  First  assume  that  Ficks  law. 


Yu  V, 
k  k 


(21) 


is  valid.  This  is  only  strictly  true  for  a  binary  mixture  where  thermal 
diffusion  is  negli^ble.  Dj^  ^  is  a  multicomponent  diffusion  coefficient. 
Equation  (21)  can  be  written  as 


2 

P 


k  m 


(22) 


As  a  first  approximation,  we  can  take 


D,  and  p  X  to  be  constant. 
K  m 

.N. 


Also,  we  can  assume  that  c  =  c  .,  k=l,2,...N.  Here  c  must  be  chosen 

p  p  K  p 

so  as  to  obtain  the  proper  adiabatic  temperature  T„.  The  procedure  for 
choosing  these  constants  and  the  justification  of  ^is  approximation  is 
given  in  reference  6.  This  level  of  approximation  is  reasonably  accurate, 
if  the  constants  are  properly  chosen. 


Equations  (15)  and  (16)  now  become 


8V, 


3t 


'J'a 


3Y. 


d>l) 


-Li 

% 


p^a 


3Y 


km 


3^ 


t„  \M^/  P  (23) 


and 


f  06  m  3T  Lss.  f  — .RAi 


'Sc^k  .^k  ,  (24) 

P  Tlo 


Equation  (5)  now  simplifies  to 


\  •  ''k°  • 

so  the  required  input  parameters  for  this  level  of  approximation  are  the 
constants  p^Dj^  p  X,  c^,  hj^°  and  the  functions  Rj^(T  ,  Yj  ,  ...»  Yj^) . 
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III.  THE  NUMERICAL  METHOD  -  PDECOL 


7 

The  package  PDECOL,  developed  by  Madsen  and  Sincovec  ,  was  used  to 
solve  the  equations.  This  package  is  designed  to  solve  a  general  system 
of  N  nonlinear  partial  differential  equations  of  at  most  second  order  on 
a  finite  interval.  In  our  coordinate  system  the  appropriate  form  is 


3u  c 

^  =  f  Ct  .  K.  .  u 


. 


(26) 


where 


-► 

u  = 

(Yj  ,  ..., 

(^ 

*■  3,|,  '  *  •  *  ’ 

”^N-1  8T  , 
3t  '  3t  ■* 

* 

(-* 

V  2  *  *  * 

3t 

2 

3  T  . 

~T  ^  • 
3t 

*  $  2  * 

3t 

(27) 


Fairly  general  boundary  conditions  of  the  form 


S  (  u,  u^)  *  \  (t) 


(28) 


are  allowed,  where  S  and  1  are  arbitrary  vector  valued  functions  with 
N  components.  Each  solution  component  is  assumed  to  be  a  known  function 
oft  at  the  initial  time  t  =  t^.  That  is,  t)»  k=l,2, . .  .N-1,  and 

T  (t  ,  t)  are  known  functions.  The  initial  conditions  must  be  consistent 
with^the  boundary  conditions. 

The  spatial  discretization  is  accomplished  by  finite  element 
collocation  methods  based  on  B-splines®.  The  user  must  supply  a  set  of 


^iV.  K.  Madsen  and  R.  F.  Sinaovea,  ''PDBOOL:  General  Collocation  SoftvKXPe 
for  Partial  Differential  Equations"^  Preprint  UCRL-78263  (Rev  I), 
Lawrence  Livermore  Laboratory ^  (1977). 

^C.  de  Boor^  "Package  for  Calculating  with  B-Splinee",  Siam,  J.  Burner. 
Anal.  ^  441-472,  (1977). 


NB  breakpoints,  that  is,  a  set  of  strictly  increasing  locations  where 
the  polynomials  are  joined.  He  must  also  supply  the  order  KORO  of  the 
splines  and  the  number  of  continuity  conditions,  NCC,  to  be  applied  at 
the  breakpoints.  That  is,  if  NCC=1,  the  approximating  function  is 
continuous;  if  NCC=2,  it  is  continuous  and  smooth;  and  so  on.  PDECOL 
then  generates  a  set  of  NC  =  KORD(NB-l)  -  NCC(NB-2)  basis  functions  and 
collocation  points.  The  basis  functions  B^(i|j)  are  piecewise  polynomials 

of  order  KORD  -  1.  The  basic  assumption  is  that  the  solution  can  be 
written  in  the  form 

N 

\  =  ^  ‘  k=l,...N  ,  (29) 

i=l 


where  the  basis  functions  B^(^)  span  the  solution  space  for  any  fixed  t 

to  within  a  small  error  tolerance.  The  time  dependent  coefficients 
c.  (i)  are  determined  uniquely  by  requiring  that  the  expansions  above 
satisfy  the  given  boundary  conditions  and  that  they  satisfy  the  partial 
differential  equations  exactly  at  the  (N-2)  interior  (collocation) 
points.  Since  by  definition  a  B-spline  is  zero  except  over  a  small 
interval,  at  any  collocation  point  no  more  than  KORD  of  the  B~splines 
are  non-zero.  So  the  system  of  ODE'S  for  the  coefficients  will 

not  be  fully  coupled. 


The  boundary  conditions  (28)  must  also  be  changed  into  ordinary 
differential  equations.  This  can  be  done  by  taking  the  derivative  with 
respect  to  time,  which  results  in 


N 

I 


9b,  9u. 

r  -I  __1 
^  9u.  9t 
1 


+ 


k»l,2...N. 


(30) 


When  eauation  (29)  is  substituted  into  equation  (30),  a  set  of  ODE's  in 
the  results. 


Unlike  a  finite  difference  code,  the  program  can  be  run  with  no 
boundary  condition  at  either  the  left  or  the  right  boundary.  The  program 
simply  collocates  at  the  boundary,  using  the  same  procedure  as  for  the 
interior  points.  It  must  be,  of  course,  the  user’s  responsibility  to 
define  a  mathematically  meaningful  PDE  problem. 
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This  system  of  ODE's  is  integrated  in  time,  using  a  variant  of  the 
Gear  stiff  integrator^.  This  is  a  fully  implicit,  pyedictor-corrector 
method.  The  required  banded  Jacobian  is  generated  internally  by  the 
program.  Once  the  integrator  has  reached  a  desired  output  time,  the 
values  of  Yj^  and  T  can  be  obtained  for  any  ij/  by  substituting  into  the 

expansions  (29) .  The  accuracy  of  the  time  integration  is  determined  by 
a  user  supplied  error  tolerance  e. 

For  the  actual  integration,  the  code  requires  the  values  of  the 
basis  functions  and  the  first  two  space  derivatives  only  at  the 
collocation  points.  Since  computing  these  values  is  complicated,  this 
is  done  once  at  the  beginning  of  the  program  and  the  values  are  saved. 

In  general,  a  user  must  supply  a  main  program  and  three  subroutines. 
The  program  MAIN  sets  the  values  of  the  required  parameters,  calls  the 
integrator,  and  writes  any  desired  output.  The  subroutine  UINIT  gives 
the  initial  conditions.^  Given  a  collocation  point  i|),  the  subroutine 

fust  return  the  vector  u  (t  ,  tji) .  The  subroutine  F  evaluates  the  function 
given  by  Eq.  (26),  Given°t,  u,  u^,  and  u^^,  the  routine  must  return 

the  vector  3u/3t.  The  program  automatically  converts  this  to  the 
corresponding  set  of  ODE's  involving  the  The  subroutine  BNDRY 

specifies  the  boundary  conditions.  That  is,  given  t,  ip,  u,  and  u.  , 
the  subroutine  returns  the  quantities  abj^/3Uj,  3bj^/3  u^  3Zj^/3t  ; 

j,k  =  1,2...  N,  for  both  and 

IV,  MODIFICATIONS  OF  PDECOL 


The  time  integration  is  controlled  by  a  user  supplied  error 
tolerance  e.  Single  step  error  estimates  divided  by  CMAX.  Ci)  will  be 

^  Ci) 

kept  less  than  e  in  the  root-mean-square  norm.  In  PDECOL,  CMAXj^ 

is  initially  set  to  the  maximum  of  Ic,  and  1.0.  Thereafter, 

CMAXj^^^^  is  the  largest  value  of 

CMAX,  if  that  is  larger.  However,  this  error  criterion  does  not 
produce  the  desired  accuracy  for  flame  simulations,  because  radical 
species  with  small  concentrations  will  control  the  flame,  and  must  be 
computed  accurately.  But  mass  fractions  much  smaller  than  1.0  will  not 
be  computed  accurately  using  the  original  PDECOL  criterion.  An  alternate 
criterion  is  the  purely  relative  error  criterion,  that  is,  CMAXj^Ci)  s 


seen  so  far,  or  the  initial 


(i) 


This  criterion  is  also  unacceptable  because  some  species 


at  some  locations  will  approach  zero,  and  excess  computation  will  result 


^A.  C.  Hindmarehf  "Prelimiruxpy  Documentation  of  GEARIB;  Solution  of 
Implicit  Syateme  of  Ordinary  Differential  Equations  with  Banded 
Jacobian" t  Rep-  UCID-30220,  Lawrence  Livermore  Laboratory t  (1976). 
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in  accurately  computing  these  negligible  concei)trations.  Consequently 
a  semi -relative  error  control  is  used.  CMAX,  is  chosen  as  the 
Cil  ^ 

maximum  of  Cj^'-  •'  and  a  user  supplied  parameter  SREC.  Thus,  mass  fractions 

less  than  SREC  will  be  computed  less  accurately.  We  have  normally  used 
SREC=10"®. 


The  original  program  PDECOL  is  fully  implicit.  That  is,  it  generates 
a  set  of  NO  =  N  X  NC  ordinary  differential  equations  of  the  form 


A  ^  =  g  (t,  c)  .  C31) 


where  c  =  .  The  resulting  Jacobian 

3g/3c  is  a  banded  NO  by  NO  matrix.  The  band  width  is  3  X  ML+1,  where 
ML  =  N  (KORD  -  13  -  1.  To  advance  a  time  step,  the  program  first 
computes  an  explicit  predictor  for  the  values  of  the  ab  the  next 

time  step.  Then  it  solves  a  set  of  linear  algebraic  equations  involving 
the  Jacobian  to  correct  the  values.  For  a  stiff  system,  such  as  one 
involving  chemistry,  this  allows  time  steps  orders  of  magnitude  larger 
than  those  of  an  explicit  method®.  The  drawback  is  the  storage  and 
execution  time  required  to  work  with  the  Jacobian. 

For  example,  consider  a  system  where  KORD  =  4,  N  =  10,  and  NC  =  20. 
The  bandwidth  of  the  Jacobian  is  88,  and  we  have  a  system  of  300  ODE's. 

We  then  have  26,400  possible  non-zero  elements  in  the  banded  matrix. 

The  amount  of  storage  required  also  increases  rapidly.  For  instance, 
if  N  3  20  instead  of  10  the  bandwidth  will  be  178,  and  we  will  require 
106,800  words  to  store  the  Jacobian.  Since  we  want  to  be  able  to  solve 
systems  at  least  this  large,  the  storage  requirements  become  almost 
prohibitive.  In  addition,  solving  such  large  systems  of  linear  algebraic 
equations  is  very  time  consuming. 

To  avoid  this  problem,  we  essentially  uncouple  the  partial  differ¬ 
ential  equations  and  solve  them  successively.  The  basic  procedure  is 
illustrated  below.  Suppose  we  are  at  time  t  and  we  want  to  advance 
a  time  step  to  We  first  integrate  the"equation  for  Yj^,  under  the 

assumption  that  Y2,...,  Yj^_j^,  T  are  constant  at  their  t^  values.  This 

uncouples  the  first  PDE  from  the  system.  Subsequently,  we  integrate  the 
second  PDE  for  Y2,  using  the  new  value  of  Yj  at  and  the  old  values 

for  the  other  variables.  Continuing  this  process,  we  finally  solve  for 
T(tn^l),  using  updated  values  for  all  the  mass  fractions. 


\ 


In  general,  this  method  is  restricted  to  smaller  step  sizes  than  a 
fully  implicit  method.  However,  as  steady  state  is  approached,  all  the 
variables  approach  constants  with  respect  to  time,  and  the  temporal 
coupling  vanishes. 

Now  consider  how  this  assumption  affects  the  associated  system  (31) 
of  ODE'S.  In  the  Jacobian,  we  will  have 

-  0  if  k  A. 

That  is,  changes  in  the  PDE  will  not  affect  the  PDE  if  k  i* 

So' most  of  the  terms  in  the  Jacobian  will  become  zero. 

The  remaining  nonzero  elements  are  not  in^banded  form.  However,  we 
can  accomplish  this  by  rearranging  the  vector  c  =  Cj(2),  ... 

^^(NC)^  That  is,  the  coefficients  for  each  PDE 

are  grouped  together  instead  of  the  coefficients  for  each  collocation 
point. 

The  Jacobian  matrix  can  now  be  decomposed  into  N  smaller  NC  X  NC 
matrices  on  the  main  diagonal.  Moreover,  the  band  width  of  these 
matrices  is  only  3XNL-»'l,  where  ML  ■  KORD  -  2.  The  Jacobian  is  essentially 
N  separate  Jacobians,  each  for  one  PDE  with  NC  collocation  points. 

The  savings  in  storage  space  is  dramatic.  For  our  previous  example 
of  KORD  ■  4,  N  s  10,  and  NC  »  30,  we  have  2100  nonzero  elements  instead 
of  26,400.  If  N  «  20,  we  have  only  4200  nonzero  elements  instead  of 
106,800. 

The  following  procedure  is  used  to  actually  integrate  a  time  step. 

The  .predictor  of  the  predictor-corrector  method  is  USPd  to  obtain  first 
estimates  for  all  the  •  Then  the  corrected  c^^i)  are  computed, 

using  the  first  small  Jacobian  and  assuming  that  the  other  do  not 

change.  Then  the  corrected  values  for  and  the  predicted  values 

for  k  >  2,  are  used,  and  the  corrected  are  computed.  This 

process  continues  through  the  N  PDE's.  The  estimated  error  is  calculated. 
If  necessary,  the  above  procedure  is  iterated.  It  is  more  efficient  to 
use  the  predicted  values  rather  than  the  values  at  the  previous  time 
step.  In  using  the  procedure,  it  is  more  efficient  to  integrate  the 
minor  species  first,  since  they  change  most  rapidly,  then  the  major 
species,  and  finally  the  temperature. 

The  above  procedure  is  similar  to  one  developed  by  Spalding  and 
Stephenson  for  use  in  a  finite  difference  code^^. 

Spalding  and  P.L.  Stephenson,  "Lcanimr  Flame  Propagation  in  Bydrogen 
+  Bromine  Mixturee",  Proo.  /?.  Soo,  Land.  A,  324,  313-337  (1971), 
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PDECOL  was  rewritten  to  use  the  successive  calculation  method.  This 
also  means  that  we  must  change  the  user  supplied  routines  F  and  BNDRY. 

Because  of  the  rearrangement  of  the  vector  c,  the  core  integrater  will 
call  F  first  at  each  collocation  point,  asking  for  the  values  of  dY^/at. 
It  will  then  repeat  for  each  PDE.  The  routine  F  is  written  so  as  to 
return  only  the  desired  time  derivative. 

Since  the  problem  can  no  longer  directly  handle  coupling  terms, 
the  boundary  conditions  (28)  must  be  uncoupled,  that  is,  the  boundary 
conditions  must  be  of  the  form 


The  rewritten  subroutine  BNDRY  evaluates  the  quantities  abj^/auj^, 

^  "  1,2...N. 

Comparisons  of  the  execution  time  of  the  fully  implicit  method  versus 
the  successive  calculation  method  have  not  been  made.  However,  our  main 
purpose  was  not  to  reduce  the  execution  time  but  to  reduce  the  storage 
requirements.  This  iterative  procedure  accomplishes  this  goal  and 
simultaneously  gives  accurate  results  in  reasonable  run  times. 

V.  ADDITIONS  TO  PDECOL 

The  computer  program  PDECOL  requires  a  set  of  user  supplied  subroutines. 
We  have  written  a  set  of  subroutines  that  casts  the  equations  into  a 
computationally  efficient  form  and  which  generates  the  required  output. 

These  routines  are  the  first  seven  listed  in  the  appendix;  namely  MAIN, 

F,  UINIT,  FLSP,  BNDRY,  BKPT,  and  RT.  In  addition,  several  auxiliary 
programs  are  mentioned,  but  without  providing  a  listing.  These  routines 
generate  sets  of  input  data  or  actual  subroutines  that  are  used  frequently, 
or  analyze  the  output  of  the  flame  code. 

The  code  can  be  run  with  several  different  options.  We  first  describe 
the  case  of  an  unbounded  flame  (NBURN  =  0)  using  the  simplest  constant 
transport  algorithm  (NTRAN  =  1) .  There  is  no  information  available  about 
the  flame  speed  or  the  solution  profiles  (NSTART  =  1) .  Later  in  this 
section  we  describe  restarting  the  integration  (NSTART  =  2) ,  and  using 
a  more  complicated  transport  algorithm  (NTRAN  >  2).  Burner  stabilized 
flames  (NBURN  >=  1)  are  discussed  in  the  next  section. 

Initially,  a  chemistry  scheme  must  be  chosen.  In  the  appendix  an 
^2  “  ®2  "  ^2  used,  with  nine  chemical  species  (H,  OH,  0,  HO2, 

H2O2,  H.,,  O2,  H2O,  N2)  and  a  set  of  thirty  reactions  involving  these 

species,  each  with  a  rate  constant  of  the  form  a  T^’  exp  (c/T)  (either  b 
or  c  can  be  zero).  This  infonnation  is  used  as  input  to  an  auxiliary 
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code.  This  code  writes  the  subroutine  RT  that  computes  the  chosistry 
terms  R.  The  procedure  is  analogous  to  that  used  to  write  the 

transport  subroutines^.  This  subroutine  is  attached  to  the  flame  code 
in  the  job  stream  and  can  be  used  for  any  problem  involving  this  set  of 
kinetics . 

The  actual  subroutine  RT  has  three  main  parts.  The  rate  constants 
are  evaluated  for  the  current  temperature  and  stored  in  the  vector  RK. 
Since  the  code  uses  a  successive  calculation  method,  the  subroutine  will 
be  called  N  times  at  each  time  step  and  each  collocation  point.  The 
rate  constants  are  only  recomputed  if  the  temperature  has  been  changed. 
Otherwise  this  section  is  skipped.  Each  rate  constant  is  multiplied 
by  the  concentrations  of  the  appropriate  reactants  divided  by  p  to 
obtain  r^/p.  The  terms  Mj^/p  (stored  in  the  vector  R)  are  calculated 

by  adding  the  rates  for  the  reactions  in  which  is  a  product, 
subtracting  the  rates  for  which  Yj^  is  a  reactant,  and  multiplying  by  M|^. 

The  choice  of  a  transport  algorithm  determines  the  form  of  the 
subroutine  F.  The  version  given  in  the  appendix  is  for  constant 
transport.  Common  statements  are  used  to  make  the  appropriate  constants 
available  for  all  subroutines.  The  chemistry  terms  required  are  found 
by  calling  RT. 


The  initial  temperature  Ty  and  mass  fractions  Y^u  of  the  unbumed 


'll 

Another  auxiliary  code  determines  the 

and  the  constants 


gas  are  input  quantities. 

adiabatic  temperatures  Tg  and  mass  fractions  Yj^g, 

p^Dkju,  h^,  pX,  and  Cp  needed  by  the  constant  transport  algorithm  (see 

reference  6).  This  information  is  saved  on  a  data  file  and  attached  to 
TAPE  11  when  the  code  is  run. 


The  remaining  data  necessary  to  run  the  program  is  read  in  on  cards 
(TAPE  5) .  The  pressure  p  and  the  normalizing  constants  t^,  and  T 
are  specified.  The  option  parameters  NSTART,  NTRAN,  and  NBURN  are  chosen. 
The  numerical  parameters  required  by  the  code  are  specified;  that  is. 


\|»,  and  i)!-,  the  final  integration  time  t. 


'L  ‘ —  "^R’  — . . .  FINAL’ 

control  parameters  e  and  SREC,  and  a  set  of  breakpoints 


the  time  integration  error 


The  major  difficulty  in  efficiently  solving  the  flame  equations  is 
choosing  an  appropriate  set  of  breakpoints.  These  must  be  close  enough 
that  spatial  errors  are  minimized  and  yet  not  so  dense  that  one's 
computer  resources  are  exceeded.  The  breakpoints  should  be  densest  in 
the  flame  front,  where  the  gradients  are  very  steep. 

The  technique  in  the  code  is  to  use  a  static  mesh,  with  the  break¬ 
points  most  closely  spaced  near  the  center  of  the  interval  of  integration. 
The  flame  front  is  then  forced  to  remain  near  the  center  of  the  interval. 
This  is  done  by  adjusting  a^,  the  mass  flux  through  the  origin.  At 
steady  state,  the  mass  flux  through  the  flame  is  a  constant.  So  m^  is 
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iteratively  modified  to  match  the  steady  state  mass  flow  through  the 
flame.  This  leads  to  a  coordinate  system  in  which  the  flame  front  is 
at  rest.  The  transient  behavior  can  cause  the  flame  ‘front  to  drift  away 
from  the  center.  This  can  also  be  corrected  by  modifying  m^. 

This  method  requires  a  procedure  for  tracking  the  flame  front.  To 
do  this  the  position  of  a  specific  temperature  Tgn  is  monitored.  The 
code  attempts  to  keep  this  temperature  located  at  the  center  of  the 
interval  of  integration.  As  a  heuristic  rule  this  temperature  is  defined 
by 

Ten  “  Tu  -  C  .4  (Tg  -  T^)  .  (33) 

The  average  of  T^  and  Tg  is  not  used.  This  is  because  the  mixture  will 
not  normally  reach  the  adiabatic  temperature  at  the  end  of  the  flame  front. 
Rather,  there  is  a  radical  overshoot,  and  the  recombination  of  these 
radicals  will  very  slowly  raise  the  temperature.  T^ji  as  defined  will 
usually  be  close  to  the  center  of  the  flame  front.  The  details  of  this 
iterative  procedure  are  discussed  in  reference  11. 


So  the  breakpoints  may  be  chosen  to  be  densest  in  the  center  of  the 
interval  of  integration.  However,  it  is  tedious  to  have  to  choose  an 
entire  breakpoint  sequence  for  each  problem.  We  have  developed  a  procedure 
to  generate  an  appropriate  type  of  breakpoint  sequence  from  a  small 
number  of  parameters.  By  varying  three  parameters,  a  wide  variety  of 
breakpoint  sequences  can  be  generated. 

The  user  must  supply  HINT,  the  number  of  intervals  (MB  »  NINT+1), 

NCN,  the  number  of  intervals  of  equal  length  that  will  be  at  the  center 
of  the  interval,  and  PC,  the  ratio  between  the  longest  intervals  (on  the 
boundaries)  and  the  shortest  intervals.  Also  let  L  be  the  total  length 
of  the  interval  of  integration,  that  is,  L  >  The  program 

generates  a  set  of  intervals  whose  lengths  increase  by  a  constant  factor 
a,  where 

a  =  log"^  [2 (log  FC)/(NINT-NCN)]  .  (34) 

The  common  length  LC  of  the  NCN  center  shortest  intervals  is 

LC  »  L/[NCN  ♦  2  ^ 

The  procedure  can  best  be  seen  by  example.  Suppose  we  have  -  0, 

^  T.P.  Coffee  and  J.M.  Heinerl,  "A  Method  for  Confuting  the  Flame  Speed 
of  a  Laminar t  Premixed j  One  Dimeneional  Flame",  BRL  Technical  Report 
ARBRL-TR-02212,  January  1980. 
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\|»_  =  10,  NINT  =  12,  NCN  *  4,  and  FC  =  6.  Then  L  *  10,  a  =  1.5651,  and 
LC  >  .3155.  The  resulting  breakpoint  sequence  is  given  in  Table  1.  Note 
that  the  4  center  intervals  are  of  the  same  length,  the  length  of  the 
intervals  then  increases  by  a  factor  of  a,  and  the  two  intervals  by  the 
boundaries  are  6  times  as  long  as  the  central  intervals.  So  the  procedure 
automatically  generates  a  set  of  breakpoints  that  are  closest  together 
near  the  center,  with  the  spacing  increasing  smoothly  toward  the 
boundaries . 

Some  experimentation  is  necessary  to  choose  the  proper  value  of  the 
above  parameters.  However,  choosing  NINT  =  12,  NCN  =  4,  and  FC  between 
4  and  8  has  worked  in  most  of  the  cases  we  have  tried.  Normally  we 
experiment  using  the  constant  transport  algorithm,  and  then  use  a  more 
realistic  transport  subroutine  once  we  have  a  good  breakpoint  sequence. 

TABLE  1.  THE  SET  OF  BREAKPOINTS  GENERATED  BY  L  =  10,  NINT  *  12, 

NCN  =  4  and  FC  *  6. 


Breakpoints 


1.8929 


3.1024 


3.8752 


4.3689 


4.6844 


5.0000 


5.3155 


5.6310 


6.1247 


6.8975 


8.1069 


10.0000 


Interval  Leneths 


1.8929 


1.2094 


.7728 


.4937 


.3155 


.3155 


.3155 


.3155 


.7728 


1.2094 


1.8929 


The  breakpoint  sequence  is  written  by  the  subroutine  BKPT.  It  also 
generates  a  larger  set  of  evaluation  points  by  interpolating  between  the 
breakpoints.  Using  Eq.  (29),  the  subroutine  VALUES  can  evaluate  the  Y. 
and  T  at  this  larger  set  of  points.  This  information  is  useful  in 
generating  detailed  output,  such  as  graphs,  or  in  performing  numerical 
integrations  (see  below) . 

The  subroutine  UINIT  writes  the  initial  profiles,  where 


'I'l  =  ’I'l  +  0-24  -  <1/^) 

'i>2  ~ 


(36) 


’'kU' 


1 

P 

^kU  *  ^^kB  -  ^kU^ 

1. 

2  \*2-*J 

1 

f 

\ 

’''l  1 

'kB  ' 

’^'r. 

The  temperature  T  is  defined  similarly. 


This  particular  definition  will  give  us  T  =  Ten  at  i|»cn  *  0.5(\|(l+'I'R) 
The  choice  of  the  particular  function  that  defines  the  and  T  between 
and  'l>2  is  not  important.  We  have  used  a  straight  line  with  success, 
but  defining  a  smooth  function  is  sli'.ghtly  more  efficient. 


A  requirement  of  PDECOL  is  that  the  initial  profiles  satisfy  the 
boundary  conditions,  in  our  case  given  by  Eqs.  (17)  and  (18).  The  above 
initial  profiles  have  the  proper  unbumed  values  at  and  are  constant 
(space  derivative  zero)  near 


For  this  case  the  subroutine  BNDRY  has  a  simple  form.  At  the 
subroutine  returns  the  values  3bj^/3u^  =  1,  3bj^/3Uj^.  ■  0,  and  3Zj^  /3t  ■  0 


At  the  conditions  are  3bj^/3u^  «  0,  3t>j^/3Uj^^l^  «  1,  and  3Zj^/3t 


0. 
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To  begin  the  integration,  a  starting  value  for  m  ,  the  mass  flux 

through  the  origin  is  required.  To  do  this,  the  code°ignores  the  time 

dependent  terms,  assumes  that  the  mass  flux  pv  is  constant,  and  uses 

Eq.  (20)  to  obtain  a  value  of  pv.  This  gives  a  reasonable  starting  value 

for  m  . 
o 

The  evaluation  of  Eq.  (20)  is  carried  out  in  the  subroutine  FLSP. 

The  integral  is  approximated  using  the  trapezoidal  rule,  where  the 
integrand  Rj^  Mj^/p  is  calculated  at  the  evaluation  points.  The  flame 

speed  Vj^(-“)  is  evaluated  for  each  species  k  and  for  several  intervals 

(a,b),  where  a  =  The  initial  value  for  is  based  on  species  N-1 

integrated  over  the  entire  interval  (i|)j^,  ij(j^) . 

The  integration  is  performed  over  several  intervals  (a,b)  because 
of  the  behavior  of  the  minor  species.  These  are  normally  close  to  zero 
at  reach  a  peak  in  the  flame  front,  and  are  close  to  zero  again  at 

For  these  species,  integrating  from  the  left  boundary  to  the  flame 

front  is  much  more  accurate. 

The  chemistry  terms  in  Eq.  (20)  are  found  by  calling  RT.  The 
diffusion  terms  for  this  case  are  found  using  Ficks  law,  Eq.  (22).  Since 
we  are  assuming  that  is  constant,  this  value  is  stored  in  the 

vector  R2D.  The  subroutine  FLSP  also  computes  the  x  values  from  the 
relation 


x  (t,  i|i)  * 


/ 


'll 


d  ip. 


(38) 


using  the  trapezoidal  rule,  and  computes  an  estimate  of  the  flame 
thickness. 

As  the  time  integration  proceeds,  control  returns  to  MAIN  at  a 
series  of  output  times.  The  present  spatial  location  of  the  temperature 
T^^  is  found,  and  this  is  used  to  find  the  average  value  of  the  mass  flux 

since  the  last  output  time.  The  function  m^(t)  is  redefined  at  these 

times  so  as  to  keep  the  flame  front  in  roughly  the  same  position.  The 
output  times  are  chosen  by  the  program  so  that  the  flame  front  will  not 
drift  too  far  between  evaluations.  Also  FLSP  is  called  so  the  user  can 
see  if  the  flame  speeds  computed  from  the  different  species  profiles  are 
approaching  a  common  value. 

At  the  final  time  writes  an  output  file.  It 

consists  of  all  the  evaluation  points  the  corresponding  x^,  the 
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values  of  Yj^  and  T,  plus  their  first  and  second  derivatives  with  respect 
to  ij).  This’^file  can  be  attached  to  an  output  routine.  By  also  attaching 
F  and  RT,  we  can  compute  and  print  out  any  quantity  in  the  steady  state 
solution  in  which  we  are  interested.  This  file  can  also  be  attached  to 
a  graphics  routine. 

Similarly,  the  program  MAIN  writes  a  restart  file.  This  consists 
of  the  present  location  of  T^^,  the  present  value  of  m  ,  and  the  values 

of  Yj^  and  T  at  the  collocation  points.  This  file  can  be  used  to  restart 

the  time  integration.  It  is  attached  to  TAPEl,  and  the  parameter  NSTART 
is  set  equal  to  2.  The  input  parameters  read  in  on  cards  can  now  be 
changed  if  desired.  UINIT  will  translate  these  input  values  to  center 
the  flame  front,  and  will  use  interpolation  to  find  the  appropriate 
values  of  the  starting  profiles  at  the  new  collocation  points.  The  old 
value  of  m^  is  used  to  start  the  integration. 

So  far  only  the  constant  transport  case  (NTRAN  =1)  has  been 
considered.  For  more  realistic  algorithms  F  is  written  by  an  auxiliary 
code  and  attached  to  the  flame  code  (NTRAN)  =  2) .  Several  different 
levels  of  approximation  can  be  used.” 

In  these  more  complicated  algorithms,  the  diffusion  velocities  V. 
are  coupled,  and  they  must  be  computed  simultaneously.  But  since  the^ 
code  uses  successive  calculation,  only  the  value  of  one  Vj^  is  required 

on  each  call  to  F.  Because  the  computation  is  time-consuming,  it  is 
preferable  not  to  compute  the  V.  N  times  at  each  time  step  for  each 
collocation  point.  To  economize  computer  time,  all  the  thermodynamic 
and  transport  quantities  required  at  all  the  collocation  points  for 
k  =  1  are  computed,  and  stored  in  vectors.  For  k  >  1,  we  use  the  same 
values,  even  though  some  of  the  Yj^  terms  have  changed  slightly.  Only 
the  chemistry  terms  are  reevaluated.  Since  the  chemistry  normally 
changes  much  more  rapidly  than  the  transport,  this  will  still  be  a  good 
approximation. 

In  generating  an  approximation  to  the  Jacobian  (using  finite 
differences)  it  is  necessary  to  recompute  the  transport  each  time  F  is 
called,  since  what  is  of  interest  is  the  effect  of  changes  in  Yj^  and  T  on 

the  time  derivatives.  Ignoring  the  changes  in  transport  leads  to  an 
inaccurate  Jacobian.  However,  the  rate  constants  k^  and  the  thermodynamic 

quantities  c^  and  hj^  need  be  recomputed  only  if  T  is  changed,  since  they 

only  depend  on  temperature. 

As  will  be  seen  in  the  next  section  it  will  still  be  useful  to 
formulate  the  mass  flux  in  a  Ficks  law  form  Eq.  (22).  For  the  more 
complicated  transport  subroutines,  we  define 


2 

This  new  quantity  P  is  no  longer  constant  with  respect  to  space  or 
time. 


VI.  BURNER  STABILIZED  FLAMES 

So  far  we  have  only  discussed  flames  that  propagate  in  an  unbounded 
medium.  In  actual  experiments,  the  flame  will  usually  be  stabilized  by 
a  burner.  It  is  useful  to  be  able  to  model  this  type  of  experiment  so 
that  we  can  compare  experimental  and  calculated  profiles. 

Our  basic  idealization  is  of  a  cylindrical,  porous  plug  burner. 

The  premixed  gas  exits  the  plug  with  a  constant  velocity  v,  but  the  plug 
prevents  back  diffusion  of  the  products  into  the  burner.  When  the  gas 
is  ignited,  a  flame  front  develops  and  propagates  toward  the  burner. 

The  gas  velocity  v  at  the  burner  must  be  less  than  the  flame  velocity, 
or  the  flame  will  be  blown  off  the  burner.  In  the  model  this  will  look 
like  an  unbounded  flame.  In  an  experiment  the  flame  will  be  extinguished 
by  the  surrounding  atmosphere.  As  the  flame  approaches  the  burner,  the 
burner  surface  acts  as  a  heat  sink  for  the  flame.  The  loss  of  heat  slows 
down  the  flame  velocity,  until  the  flame  stabilizes  near  the  surface  of 
the  burner  at  the  gas  velocity.  As  the  gas  velocity  v  is  decreased,  the 
flame  loses  more  heat  to  the  burner,  and  stabilizes  closer  to  the  burner 
surface.  If  v  is  made  too  small  in  an  experiment,  the  flame  can  flash 
back  into  the  burner. 


For  this  kind  of  burner,  air  is  entrained  along  the  outside  edges 
of  the  flame.  However,  the  center  of  the  flame  will  correspond  closely 
to  a  premixed,  laminar,  one-dimensional  flame. 

The  only  change  in  the  equations  is  in  the  left  boundary  condition 
at  the  surface  of  the  burner.  Because  of  back  diffusion  from  the  flame 
to  the  burner  surface,  the  mass  fractions  of  the  unburned  mixture  are 
not  conserved.  However,  due  to  conservation  of  mass,  the  mass  flux 
fractions,  defined  as 


,  .  \  \  ''k 

k  pv  ”^k*pv 


(40) 


are  conserved. 

‘kU  ■  ''kU- 
surface  are 


Within  the  burner,  the  diffusion  velocities  are  zero. 
So  the  appropriate  boundary  conditions  at  the  burner 


E 


k 


kU 


(41) 


The  boundary  conditions  for  the  temperature  equation  depends  on  how  heat 
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is  extracted.  We  assume  that  the  burner  is  maintained  at  a  constant 
temperature.  Then  the  boundary  condition  is  the  same  as  for  an  unbounded 
flame. 


T  =  Ty.  (42) 

12 

This  idealization  is  discussed  by  Hirshfelder,  Curtiss  and  Bird 

To  implement  this  in  the  flame  code,  the  input  parameter  NBURN 
is  set  equal  to  1.  The  fluid  velocity  at  the  burner  must  be  specified. 
Then  m  is  a  predetermined  constant  instead  of  an  adjustable  parameter. 
All  thS  other  input  data  remains  the  same.  However,  the  code  will  handle 
this  data  differently. 

The  breakpoints  are  generated  differently,  since  the  flame  front 
will  be  near  the  left  boundary  instead  of  in  the  center  of  the  interval 
of  integration.  NCN  is  the  number  of  intervals  of  equal  length  at  the 
left  boundary.  FC  is  the  ratio  between  the  longest  interval  (at  the 
right  boundary)  and  the  shortest  interval  (at  the  left  boundary) .  The 
constant  a  is  now  defined  by 

a  =  log'^  [(log  FC)/(NINT-NCN)].  (43) 


The  common  length  LC  of  the  NCN  shortest  intervals  is  now 

LC  =  L/[NCN  +  a(a*^'^^^'^^^^-l)/(a-l)].  (44) 

As  an  example,  consider  the  same  input  parameters  that  we  used  for 
an  unbounded  flame;  ij;,  =  0,  =  10,  NINT  =  12,  NCN  =4,  and  FC  =  6. 

Li  K 

Then  L  =  10,  a  =  1.2510  and  LC  =  .4181.  The  resulting  breakpoint  sequence 
is  given  in  Table  2.  The  shortest  intervals  are  near  the  burner  surface, 
where  fairly  steep  gradients  are  expected. 

In  general,  it  is  possible  to  choose  a  shorter  interval  of  integration 
L  for  a  burner  stabilized  flame.  There  is  no  longer  a  need  for  a 
relatively  long  interval  to  the  left  of  the  flame  front  in  order  to 
approach  the  unburned  conditions.  The  number  of  breakpoints  can  then 
also  be  reduced. 

The  initial  profiles  are  also  chosen  differently.  The  code  sets 
^J.O.  Hirshfelder,  C.F.  Curtiss,  and  R.B.  Bird,  op.  ait.,  pp.  761-763. 


TABLE  2.  THE  SET  OF  BREAKPOINTS  GENERATED  BY  L-10,  NINT  »  12, 
NCN  =  4  and  FC  =  6  (NBURN  =  11 

Breakpoints  Interval  Lengths 


.3458 

.3458 

.3458 

.6916 

.3458 

1.0374 

.3458 

1.3832 

.4326 

1.8159 

.5412 

2.3571 

.6771 

3.0342 

.8471 

3.8812 

1.0597 

4.9409 

1.3257 

6.2666 

1.6585 

7.9251 

2.0749 

10.0000 

h  ^^L  ^  0.1 

'‘'2  =  \  ^  0.3 
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and  then  uses  Eq.  (37)  as  before.  This  flame  front  is  more  likely  to  be 
close  to  the  position  of  the  final  steady  state  stabilized  flame.  Also, 
the  flame  front  is  chosen  to  be  narrower.  If  Eq.  (36)  is  used,  the 
initial  transient  flame  velocity  may  be  less  than  the  velocity  of  the 
fluid.  The  flame  then  drifts  toward  the  right,  and  can  be  carried 
completely  outside  the  interval  of  integration  before  it  stabilizes  . 
Choosing  a  steeper  flame  front  leads  to  a  larger  initial  flame  velocity. 

The  boundary  condition,  Eq.  (41),  can  cause  some  difficulty.  It  is 
fairly  straightforward  for  the  constant  transport  case.  Using  Fick's 
law,  Eq.  (22),  it  can  be  rewritten  as 
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Recall  that  PDECOL  converts  this  to  a  time  dependent  equation  hy  taking 
the  time  derivative.  In  the  form  used  by  PDECOL,  Eq.  (46)  becomes 


So  the  subroutine  BNDRY  must  now  return  the  values  3bj^/3Uj^  =  1  , 

3bj^/3Uj^  a  -  ™o’  =  0  at  for  k  =  1,2...  N  -  1. 

The  mass  flux  fractions  will  then  have  the  proper  values  as  the 
integration  proceeds. 

For  more  realistic  transport  algorithms  the  mass  flux  p  Yj^  is  a 

complicated2 function  of  all  the  mass  fractions.  But  by  using  Eq.  (39) 
to  define  Dj^  ,  the  boundary  condition  can  be  put  in  the  same  form  as 
Eq.  (46),  excep?  that p ^D.  is  not  a  constant.  Taking  the  time  derivative 
we  obtain 


L.  r  ^ 

3t  '•  3ij<  ^  " 


I  »  0.  (48) 
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The  time  derivative  of  p  Dj^  can  not  be  evaluated.  It  is  not 

possible  to  write  it  as  an  explicit  function  of  the  Yj^  and 

Moreover,  since  the  code  uses  successive  calculation,  cross  coupling 
terms  between  the  different  species  are  not  allowed.  By  necessity,  the 
last  term  in  Eq.  (48)  must  be  ignored. 


It  is  possible  to  just  use  the  expression 


3t 


(49) 


as  the  boundary  condition.  At  steady  state,  the  mass  flux  fractions2 
will  approach  constants.  But  because  of  the  transient  behavior  of  p 
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the  mass  flux  fractions  will  converge  to  incorrect  values. 


What  is  required  is  a  correction  term  that  will  cause  the  code  to 
converge  to  the  proper  steady  state  limit.  Moreover,  it  must  be  of  a 
form  that  can  be  used  in  the  code.  This  can  be  done  by  choosing  the 
boundary  condition 


”'k  “km  a  § 

at  ■  I|/^m^  at  ^a.p  ’^l 


(50) 


where 


at 


^kU  ■ 

Li 

^FINAL 


(51) 


The  rationale  is  that  changes  rather  slowly.  As  it  changes  the 

mass  flux  fraction  at  <11^^  will  change  slowly  from  the  proper  value.  We 

use  a  heuristically  defined  function  az/at  to  modify  the  value  of  e. 
until  it  again  equals  Y.  ,  So  as  the  time  integration  proceeds,  the* 
value  of  at  will  vary  slightly,  but  will  approach  the  proper 

steady  state  limit. 

It  is  necessary  to  check  that  the  boundary  conditions  are  still 
consistent  with  the  initial  conditions.  For  the  initial  profiles  defined 
by  Eq.  (37)  this  will  be  the  case.  The  mass  fractions  Yj^  are  constant 

near  so  3Yj^/3i|)  at  will  be  zero.  Also  3T/3\1»  at  is  zero.  For 

any  transport  algorithm  all  the  will  be  zero  at  and  Cj^  •  Yj^.  So 

the  initial  profile  is  consistent  with  the  boundary  conditions. 

When  the  integration  is  restarted  (NSTART  =  2),  the  space  derivatives 
at  are  normally  non-zero.  It  is  necessary  to  modify  the  way  PDECOL 

determines  the  initial  values  of  Normally  this  is  done  by  the 

subroutine  INITAL.  It  calls  the  user  supplied  subroutine  UINIT  to 
obtain  the  values  Uj^Ct^  ,  ) »  where  the  j  •  1,2...  NC  are  the 

collocation  points.  Then  by  substituting  into  the  expansion  (29),  it 
obtains  a  set  of  linear  algebraic  equations  of  the  form 

NC 

u,^(to.  'I'j)  *  I  ‘'k^^CV  ^  *  1,2. ..N  .  j  -  1,2. ..NC  .  (52) 

i*l 
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These  systems  of  equations  are  solved  to  obtain  the  initial  Cj^ 


(i) 


For  a  burner  stabilized  flame,  the  appropriate  vhlues  of  the  mass 
fractions  at  i|;i  =  41^  are  not  known  but  the  values  of  the  mass  flux 
fractions  e,  are.  Thus,  we  use  the  boundary  condition  (46),  substitute 
the  expansion  (29) ,  and  obtain 


Y 


kU 


Bi(tPi) 


3B.(4/j) 

Hjj 


k=l,2...N-l.  (53) 
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The  values  of  p  Dj^  from  the  last  nan  of  the  code  are  used.  These 
values  are  also  sai^d  in  the  restart  file.  So  for  a  burner  stabilized 
flame,  the  subroutine  INITAL  has  been  changed  so  that  it  generates  this 
new  set  of  equations  for  j  =  1,  k  =  1,2...N  -  1. 

Note  that  the  above  procedure  is  not  strictly  necessary.  Since  we 
have  a  correction  teim  3Z/3t,  the  boundary  conditions  will  approach  the 
proper  values,  even  if  they  are  incorrect  at  the  start  of  the  integration. 
But  it  is  more  efficient  to  begin  with  the  boundary  conditions  as  accurate 
as  possible. 

To  actually  make  comparisons  with  experiments,  some  of  the  boundary 
conditions  may  have  to  be  further  modified.  For  instance,  evidence  exists 
that  hydrogen  atoms  H  combine  very  rapidly  on  the  buimer  surface  to  form 
H2.^^  For  practical  purposes,  we  consider  this  to  be  instantaneous.  Then 
if  our  first  species  is  H  and  our  second  species  is  H2,  their  boundary 
conditions  are 


Yj  =  0  at  iI'l  C54) 

and' 


=1  ^  *=2  “  ^2u  ’•'l  ' 

These  changes  can  easily  be  made  in  UINIT  and  INITAL. 

As  the  time  integration  proceeds,  control  returns  to  the  routine 
MAIN  at  five  equally  spaced  output  times.  The  code  calls  subroutine 
FLSP  so  the  integration  process  can  be  monitored.  There  is  no  need  to 
adjust  m  ,  since  this  is  now  a  predetermined  constant. 
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30 


VII.  NIMERICAL  CCMSIDERATIONS 


This  code  has  been  applied  to  the  H2-O2-N.  system,  with  nine  chemical 
species.”  An  earlier  version  of  the  code  has  Seen  applied  to  the  ozone 
system,  with  three  species^^.  Several  calculations  have  been  performed 
for  methane-air  flames.  In  each  case  the  code  integrates  in  time  until 
the  steady  state  solution  has  been  reached.  What  is  desired  is  accurate 
values  for  the  flame  speed  and  accurate  temperature  and  species  profiles. 
Results  are  not  discussed  here,  since  they  are  given  in  the  papers 
referenced  above.  Instead,  we  will  discuss  some  of  the  considerations 
necessary  to  obtain  accurate  results,  based  on  our  experiences  with  the 
above  systems. 

The  flame  speed  can  be  calculated  from  any  of  the  species  profiles 
using  Eq.  (20).  Normally  the  time  integration  proceeds  until  the  flame 
speeds  based  on  the  major  species  (reactants  and  products)  agree  to 
within  a  small  fraction  of  a  percent.  At  this  point  the  solution  is 
very  close  to  steady  state.  The  values  of  the  flame  speed  calculated 
from  the  minor  species  (radicals)  are  very  sensitive,  and  require  more 
accuracy  (both  spatial  and  temporal)  in  order  to  achieve  the  same 
agreement . 

In  general,  it  is  easier  to  obtain  an  accurate  solution  for  a  fast 
flame  than  for  a  slow  flame.  The  slow  flame  requires  many  more 
integration  steps  before  all  the  oscillations  die  out  and  a  steady  state 
solution  is  achieved. 

To  obtain  an  accurate  solution,  both  the  spatial  and  temporal  accuracy 
must  be  sufficient.  If  the  temporal  accuracy  is  low,  the  calculated  flame 
speeds  can  oscillate  around  the  correct  value.  If  the  spatial  accuracy 
is  low,  the  problem  may  come  very  close  to  convergence,  but  to  the  wrong 
value,  since  the  communication  between  the  collocation  points  is 
inadequate.  As  a  consequence  of  the  spatial  or  temporal  accuracy  becoming 
too  small,  the  integration  may  break  down  completely. 

Spatial  accuracy  is  more  important  in  the  flame  front,  where  we 
have  steep  gradients.  Hence,  our  algorithms  for  choosing  breakpoints 
concentrates  them  in  the  region  to  be  occupied  by  the  flame  front. 

The  need  for  spatial  and  temporal  accuracy  is  connected.  Suppose 
the  number  of  breakpoints  is  increased,  but  the  temporal  error  tolerance 
e  is  not  decreased.  The  calculated  solutions  will  oscillate,  and  eventually 
the  solution  will  break  up.  On  the  other  hand,  suppose  c  is  decreased 
but  the  breakpoint  sequence  is  not  changed.  We  can  still  obtain  an 
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answer,  but  the  integration  will  take  longer,  and  the  accuracy  of  the 
solution  is  not  increased.  At  this  stage,  the  appropriate  value  of  e 
for  a  given  breakpoint  sequence  is  a  matter  of  trial  and  error. 

An  additional  problem  can  occur  with  very  slow  flames.  If  the 
interval  of  integration  is  chosen  too  small,  there  may  be  a  noticeable 
gradient  in  the  temperature  profile  at  the  cold  boundary.  Since  the  gas 
is  entering  through  the  cold  boundary  quite  slowly,  this  can  result  in 
a  substantial  heat  loss  through  the  boundary.  This  heat  loss  will  slow 
down  the  flame.  The  integration  will  converge,  but  the  computed  flame 
speed  will  be  too  low.  This  effect  can  occur  for  any  flame.  However, 
for  most  flames,  a  noticeable  heat  loss  is  due  to  a  fairly  steep  tempera¬ 
ture  gradient  at  the  boundary.  It  is  then  obvious  that  the  interval 
of  integration  must  be  increased.  But  for  a  slow  flame  (less  than  20 
cm/sec)  a  noticeable  error  in  the  flame  speed  can  result  (5  to  10  percent) , 
even  when  the  gradient  at  the  cold  boundary  appears  negligible. 

Most  of  the  cases  we  have  run  have  been  with  HINT  *  12,  KORD  «  4.  and 
NCC  =  2  (13  breakpoints,  26  collocation  points).  We  have  let  e  =  10-^, 
and  SREC  =  10"6.  The  number  of  central  intervals  NCN  was  4,  and  the 
ratio  FC  between  the  longest  and  shortest  intervals  was  between  4  and  8. 

The  length  of  the  interval  of  integration  L,  as  well  as  the  optimum 
value  for  FC,  is  a  matter  of  trial  and  error.  However,  a  reasonable  set 
of  breakpoints  can  be  found  fairly  easily  by  experimenting  with  the 
above  scheme,  using  the  constant  transport  subroutine.  The  code  can  then 
be  run  with  a  more  realistic  transport  algorithm,  using  the  restart 
option. 

To  make  sure  that  convergence  was  obtained,  the  code  was  also  run 
with  NINT  =  16  and  e  =  3  x  10-^  for  a  number  of  cases.  Only  negligible 
differences  occurred. 

The  program  was  run  on  the  BRL  CYBER  76.  For  hydrogen-oxygen 
flames,  the  run  time  varied  from  about  30  seconds  with  the  simpliest 
transport  subroutine  to  5  minutes  or  more  for  the  most  complicated  case. 

An  earlier  version  of  this  code  was  applied  to  the  solution  of  an 
unbounded  ozone  flame^^.  The  successive  calculation  procedure  was  not 
implemented  at  that  time.  However,  since  the  number  of  species  is  so 
small,  the  savings  in  reducing  the  size  of  the  Jacobian  is  not  very 
important. 
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At  that  time  we  reported  using  NINT  ■  57  to  NINT  *  70,  KORD  ■  6 
and  NCC  «  5.  The  number  of  collocation  points  was  then  between  62  and 
75.  We  have  now  been  able  to  solve  this  case  with  equal  accuracy  using 
NINT  «  12,  KORD  «  4,  and  NCC  ■  2.  Partly  this  is  due  to  choosing  NCC  •  2 
The  program  then  generates  a  special  choice  of  collocation  points 
(Gauss  -  Legendre  quadrature  points  in  each  subinterval)  which  gives 
increased  spatial  accuracy.  More  important  is  the  improved  algorithm 
for  choosing  breakpoints.  By  using  a  little  more  care,  a  much  smaller 
set  of  breakpoints  can  still  reproduce  the  flame  with  sufficient  accuracy 
So  the  careful  choice  of  breakpoints  is  probably  the  most  important 
consideration  in  running  the  code  efficiently. 
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APPENDIX  A 


APPENDIX  A 


A  listing  of  the  computer  code  follows.  The  subroutines  MAIN,  F, 
UINIT,  FLSP,  BNDRY,  BKPT,  and  RT  are  written  for  laminar  flame  problems. 
The  subroutine  F  given  is  for  the  constant  transport  assumption.  The 
subroutine  RT  given  is  for  a  H2-O2-N2  flame.  There  are  nine  chemical 

species  and  thirty  reactions. 

The  rest  of  the  subroutines  are  from  PDECOL.  They  have  been 
modified  as  discussed  in  the  text. 

After  the  listing,  the  job  stream  and  output  is  given  for  a  typical 
hydrogen-oxygen  flame.  The  initial  unburned  gas  is  50%  H-  and  50%  air, 
where  air  is  21%  O2  and  79%  N2.  The  nitrogen  is  considered  to  be  a 
diluent,  and  does  not  react.  The  required  starting  information  is 
attached  to  TAPEll.  The  main  program  and  the  subroutine  RT  are  attached 
in  a  compiled  form.  The  integration  begins  from  the  initial  profiles 
(37) .  At  the  final  output  time,  the  flame  speed  and  the  temperature  and 
species  profiles  are  given,  as  well  as  the  corresponding  x  values.  The 
restart  file  is  written  on  TAPE  2,  and  the  output  file  on  TAPE  9.  Both 
files  are  catalogued  for  possible  future  use. 
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GLC»SARY 


Cpj^  ■  specific  heat  of  species  k,  cal-g"^  -  K'^. 

Cp  ■  specific  heat  of  the  mixture,  cal-g”^  -  k"^. 

=  diffusion  coefficient  of  species  k  in  the  mixture,  cm^-s'^. 
hj^  =  specific  enthalpy  of  species  k,  cal-g"^. 

kj  =  rate  constant  for  reaction  j  in  centimeter  -  mole  -  second  units. 

-2  -1 

m^  -  mass  flux  of  the  mixture  through  the  origin,  g-cm  -s 

-  molecular  weight  of  species  k,  g-mole~^. 

N  =  number  of  chemical  species  (also  number  of  PDE's). 

NO  =  number  of  ODE's. 

p  =  pressure,  atmos. 

r^  =  rate  of  reaction  j,  mole-cm"^-s'^. 

3  -1  -1 

R  =  gas  constant  »  82.05  cm  -atoms-mole  -K  . 

Rj^  *  rate  of  production  of  species  k,  mole-cm~^-s"^. 

Sy  =  velocity  of  the  flame  relative  to  the  unburhed  mixture,  cm-s 

A  A  A 

t  3  temporal  coordinate,  s.  (t,  x  coordinate  system). 

-  temporal  coordinate,  s.  (t,  tjl  coordinate  system). 

t  =  non  dimensional  temporal  coordinate,  with  respect  to  t^  (t,  ^ 
coordinate  system) . 

A 

T  -  temperature  of  the  mixture,  K. 

^  *  temperature  of  the  mixture,  K. 

T  =  non  dimensional  temperature  of  the  mixture,  with  respect  to  T^. 

"  temperature  of  the  unbumed  mixture. 
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GLOSSARY  (continued) 

Tq  b  temperature  of  the  burned  mixture  (adiabatic  temperature) . 

D 

V  a  fluid  velocity  of  the  mixture,  cm-s~^. 

Vj^  =  diffusion  velocity  of  species  k,  cm-s"^. 

A 

X  a  Spatial  coordinate,  cm. 

Y|^  a  mass  fraction  of  species  k. 

Yj^y  a  mass  fraction  of  species  k  in  the  unbumed  mixture. 

Yy^g  a  mass  fraction  of  species  k  in  the  burned  mixture. 

a  mass  flux  fraction  of  species  k. 

X  a  heat  conductivity  of  the  mixture,  cal-cm~^-s"^-K 

.3 

p  a  density  of  the  mixture,  g-cm  . 

^  a  transformed  distance  coordinate,  g-cm'^. 

t|»  a  non  dimensional  transformed  distance  coordinate,  with  respect  to 
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